Los datos que se van a analizar en este proyecto han sido obtenidos desde Kaggle. Contienen precios de casas que fueron vendidas desde mayo de 2014 hasta mayo de 2015 en King County que es un condado ubicado en el estado estadounidense de Washington.
Lo que queremos hacer con estos datos es predecir el precio de las casas dependiendo de los datos recogidos.
En nuestro estudio inicial, la variable “Precio” es continua, por lo que vamos a categorizarla. Se van a realizar dos tipos de cetegorizaciones:
Para decidir las categorizaciones se ha optado por los cuantiles.
#Categorizamos la variable respuesta price:
quantile(datos$price, prob=seq(0, 1, length = 5))
## 0% 25% 50% 75% 100%
## 75000 321950 450000 645000 7700000
datos$price_categ1 <- cut(datos$price, breaks = c(0, 500000, 100000000), labels = c("B1", "C1"))
table(datos$price_categ1)
##
## B1 C1
## 12560 9053
datos$price_categ2 <- cut(datos$price, breaks = c(0, 330000, 650000, 100000000), labels = c("B2","M2", "C2"))
table(datos$price_categ2)
##
## B2 M2 C2
## 5881 10525 5207
A continuación se muestran cómo están categorizados los datos:
En este mapa se puede visualizar cómo se distribuyen las casas “Caras” y “Baratas”. Se observa que las casas que están cercanas al agua y a cerca de Seattle por la parte Norte, son más caras y hacia el sur más baratas.
center_lon = median(datos$long,na.rm = TRUE)
center_lat = median(datos$lat,na.rm = TRUE)
factpal2 <- colorFactor(c("green","red"),
datos$price_categ1 )
leaflet(datos) %>% addProviderTiles("Esri.NatGeoWorldMap") %>%
addCircles(lng = ~long, lat = ~lat,
color = ~factpal2(datos$price_categ1)) %>%
# controls
setView(lng=center_lon, lat=center_lat,zoom = 12) %>%
addLegend("bottomright", pal = factpal2 , values = ~datos$price_categ1,
title = "Tipos de Casas",
opacity = 1)
En este otro mapa se puede visualizar cómo se distribuyen las casas según el precio en tres categorías:“Caras”, Medio y “Baratas”.
center_lon = median(datos$long,na.rm = TRUE)
center_lat = median(datos$lat,na.rm = TRUE)
factpal2 <- colorFactor(c("green","red","yellow"),
datos$price_categ2 )
leaflet(datos) %>% addProviderTiles("Esri.NatGeoWorldMap") %>%
addCircles(lng = ~long, lat = ~lat,
color = ~factpal2(datos$price_categ2)) %>%
# controls
setView(lng=center_lon, lat=center_lat,zoom = 12) %>%
addLegend("bottomright", pal = factpal2 , values = ~datos$price_categ2,
title = "Tipos de Casas 3 categorías",
opacity = 1)
Se va a separar los datos en los 3 conjuntos de datos fundamentales:
num_total=nrow(datos)
set.seed(122556) #reproductividad
# 70% para train
indices_train = sample(1:num_total, .7*num_total)
datos_train = datos[indices_train,]
# 15% para test
indices=seq(1:num_total)
indices_test=indices[-indices_train]
indices_test1 = sample(indices_test, .15*num_total)
datos_test = datos[indices_test1,]
# 15% para validacion
indices_validacion=indices[c(-indices_train,-indices_test1)]
datos_validacion=datos[indices_validacion,]
Se van a realizar transformaciones de un conjunto de variables, estas transformaciones se aplicarán a cada conjunto de datos, train, test y validación:
Se realiza una transformación logarítmica sobre las variables sqft_living (pies cuadrados de la casa), sqft_lot (pies cuadrados del jardín) y sqft_above (pies cuadrados poe encima del suelo), hay que aclarar que esta última variable esla diferencia entre sqft_living y sqft_basement por lo que va a estar altamentente correlad con sqft_living.
Se categorizan las variables:
Bathroom, esta varible puede tomar valores decimales de 0.25 en 0.25. El número de baños se contabiliza por las piezas y cada baño completo tiene 4 piezas, por lo que con la nueva agrupación toma valores de 1 a 8 baños.
Sqft_basement, se categoriza como 0 las casas que no tienen sótano y 1 las casas que sí tienen sótano.
Grade, se va a categorizar del siguiente modo con valor 0=calidad Baja, 1= calidad media y 2= calidad alta.
Year_renovated, se categoriza como 0 = no ha tenido renovación y 1 = sí ha tenido renovación.
Se pasan a factor las variables: waterfront, view, condition, grade_categ y zipcode
Se eliminan Outliers.
datos_train <- datos_train[,-2]
datos_train$id <- as.factor(datos_train$id)
datos_train$bathrooms_group <- cut(datos_train$bathrooms,breaks = c(-1,0.25,1,2,3,4,5,6,7,8),labels=c(0,1,2,3,4,5,6,7,8))
datos_train$bathrooms_group <- as.numeric(as.character(datos_train$bathrooms_group))
datos_train$log_sqft_living <- log10(datos_train$sqft_living)
datos_train$log_lot <- log10(datos_train$sqft_lot)
datos_train$log_above <- log10(datos_train$sqft_above)
datos_train$sqft_basement_cat <- cut(datos_train$sqft_basement,breaks = c(-1,0,6000),labels=c(0,1))
datos_train$waterfront<-as.factor(datos_train$waterfront)
datos_train$view<-as.factor(datos_train$view)
datos_train$condition<-as.factor(datos_train$condition)
datos_train$grade_categ <- cut(datos_train$grade, breaks = c(0,4,9,13), labels = c(0,1,2))
datos_train$yr_renovated_catg <-cut(datos_train$yr_renovated, breaks=c(-0.5,1933, 2015), labels= c("0","1"))
datos_train$zipcode<-as.factor(datos_train$zipcode)
datos_train$posicion<-c(1:nrow(datos_train))
indices_cero_habitaciones<-datos_train[datos_train$bedrooms==0,]$posicion
datos_train<-datos_train[-indices_cero_habitaciones,]
datos_train$posicion<-c(1:nrow(datos_train))
indices_cero_banos<-datos_train$posicion[datos_train$bathrooms_group==0]
datos_train<-datos_train[-indices_cero_banos,]
datos_train$posicion<-c(1:nrow(datos_train))
indice_hab33 <- datos_train[datos_train$bedrooms==33,]$posicion
datos_train[datos_train$posicion == indice_hab33,]$bedrooms = 3
Se han realizado dos categorizaciones adicionales sobre el conjunto de datos. En este caso para ver cómo clasificaban estas variables se ha usado un arbol de decisión, las variables son: zipcode y para bathrooms_group
model_selec_zipcode<-rpart(price_categ1~zipcode,data=datos_train ,parms=list(split="gini"))
print(model_selec_zipcode)
## n= 15119
##
## node), split, n, loss, yval, (yprob)
## * denotes terminal node
##
## 1) root 15119 6385 B1 (0.5776837 0.4223163)
## 2) zipcode=98001,98002,98003,98010,98011,98014,98019,98022,98023,98024,98028,98030,98031,98032,98034,98038,98042,98045,98055,98056,98058,98059,98070,98092,98106,98108,98118,98125,98126,98133,98144,98146,98148,98155,98166,98168,98178,98188,98198 8243 1336 B1 (0.8379231 0.1620769) *
## 3) zipcode=98004,98005,98006,98007,98008,98027,98029,98033,98039,98040,98052,98053,98065,98072,98074,98075,98077,98102,98103,98105,98107,98109,98112,98115,98116,98117,98119,98122,98136,98177,98199 6876 1827 C1 (0.2657068 0.7342932) *
Se va a categorizar en dos zona1 y Zona2.
datos_train$zona<-recode(datos_train$zipcode, "98001=1; 98002=1; 98003=1; 98010=1; 98011=1; 98014=1; 98019=1; 98022=1; 98023=1; 98024=1; 98028=1; 98030=1; 98031=1; 98032=1; 98034=1; 98038=1; 98042=1; 98045=1; 98055=1; 98056=1; 98058=1; 98059=1; 98070=1; 98092=1 ;98106=1; 98108=1; 98118=1; 98125=1; 98126=1; 98133=1; 98144=1; 98146=1; 98148=1; 98155=1; 98166=1; 98168=1; 98178=1; 98188=1; 98198=1; 98004=2; 98005=2; 98006=2; 98007=2; 98008=2; 98027=2; 98029=2; 98033=2; 98039=2; 98040=2; 98052=2; 98053=2; 98065=2; 98072=2; 98074=2; 98075=2; 98077=2; 98102=2; 98103=2; 98105=2; 98107=2; 98109=2; 98112=2; 98115=2; 98116=2; 98117=2; 98119=2; 98122=2; 98136=2; 98177=2; 98199=2")
datos_train$zipcode = NULL
model_selec_bathrooms<-rpart(price_categ1~bathrooms_group,data=datos_train ,parms=list(split="gini"))
print(model_selec_bathrooms)
## n= 15119
##
## node), split, n, loss, yval, (yprob)
## * denotes terminal node
##
## 1) root 15119 6385 B1 (0.5776837 0.4223163)
## 2) bathrooms_group< 2.5 7282 1850 B1 (0.7459489 0.2540511) *
## 3) bathrooms_group>=2.5 7837 3302 C1 (0.4213347 0.5786653) *
En el resultado del modelo se ve que corta en el número de baños <2.5, por lo que se va a categorizar como 0 aquellas casas que tengan de 1 a 2 baños y como 1 las casas que tengan más de 2 baños.
datos_train$bathrooms_group <- cut(datos_train$bathrooms_group, breaks = c(-1,2.5,8),labels=c(0,1))
#0 1ó 2 baños, 1 de 3 baños en adelante
datos_train_limpio <- datos_train[c(3,22:26,8:10,27,16,17,30,28,20,21)]
#Eliminamos sqft_above porque es una combinalción lineal de sqft_living, están altamente correladas........
datos_train_limpio$log_above = NULL
datos_train_numeric <- datos_train_limpio %>% select_if(is.numeric)
Realizamos todas las transformaciones y categorizaciones para el conjunto de datos de test.
datos_test <- datos_test[,-2]
datos_test$id <- as.factor(datos_test$id)
datos_test$bathrooms_group <- cut(datos_test$bathrooms,breaks = c(-1,0.25,1,2,3,4,5,6,7,8),labels=c(0,1,2,3,4,5,6,7,8))
datos_test$bathrooms_group <- as.numeric(as.character(datos_test$bathrooms_group))
datos_test$log_sqft_living <- log10(datos_test$sqft_living)
datos_test$log_lot <- log10(datos_test$sqft_lot)
datos_test$log_above <- log10(datos_test$sqft_above)
datos_test$sqft_basement_cat <- cut(datos_test$sqft_basement,breaks = c(-1,0,6000),labels=c(0,1))
datos_test$waterfront<-as.factor(datos_test$waterfront)
datos_test$view<-as.factor(datos_test$view)
datos_test$condition<-as.factor(datos_test$condition)
datos_test$grade_categ <- cut(datos_test$grade, breaks = c(0,4,9,13), labels = c(0,1,2))
datos_test$yr_renovated_catg <-cut(datos_test$yr_renovated, breaks=c(-0.5,1933, 2015), labels= c("0","1"))
datos_test$zipcode<-as.factor(datos_test$zipcode)
#codificar la variable Zipcode
datos_test$zona<-recode(datos_test$zipcode, " 98001=1; 98002=1; 98003=1; 98010=1; 98011=1; 98014=1; 98019=1; 98022=1; 98023=1; 98024=1; 98028=1; 98030=1; 98031=1; 98032=1; 98034=1; 98038=1; 98042=1; 98045=1; 98055=1; 98056=1; 98058=1; 98059=1; 98070=1; 98092=1 ;98106=1; 98108=1; 98118=1; 98125=1; 98126=1; 98133=1; 98144=1; 98146=1; 98148=1; 98155=1; 98166=1; 98168=1; 98178=1; 98188=1; 98198=1; 98004=2; 98005=2; 98006=2; 98007=2; 98008=2; 98027=2; 98029=2; 98033=2; 98039=2; 98040=2; 98052=2; 98053=2; 98065=2; 98072=2; 98074=2; 98075=2; 98077=2; 98102=2; 98103=2; 98105=2; 98107=2; 98109=2; 98112=2; 98115=2; 98116=2; 98117=2; 98119=2; 98122=2; 98136=2; 98177=2; 98199=2")
datos_test$zipcode = NULL
datos_test$bathrooms_group <- cut(datos_test$bathrooms_group, breaks = c(-1,2.5,8),labels=c(0,1))
#0 1ó 2 baños, 1 de 3 baños en adelante
datos_test_limpio <- datos_test[c(3,22:26,8:10,27,16,17,29,28,20,21)]
datos_test_limpio$log_above = NULL
datos_test_numeric <- datos_test_limpio %>% select_if(is.numeric)
Realizamos todas las transformaciones y categorizaciones para el conjunto de datos de Validación.
datos_validacion <- datos_validacion[,-2]
datos_validacion$id <- as.factor(datos_validacion$id)
datos_validacion$bathrooms_group <- cut(datos_validacion$bathrooms,breaks = c(-1,0.25,1,2,3,4,5,6,7,8),labels=c(0,1,2,3,4,5,6,7,8))
datos_validacion$bathrooms_group <- as.numeric(as.character(datos_validacion$bathrooms_group))
datos_validacion$log_sqft_living <- log10(datos_validacion$sqft_living)
datos_validacion$log_lot <- log10(datos_validacion$sqft_lot)
datos_validacion$log_above <- log10(datos_validacion$sqft_above)
datos_validacion$sqft_basement_cat <- cut(datos_validacion$sqft_basement,breaks = c(-1,0,6000),labels=c(0,1))
datos_validacion$waterfront<-as.factor(datos_validacion$waterfront)
datos_validacion$view<-as.factor(datos_validacion$view)
datos_validacion$condition<-as.factor(datos_validacion$condition)
datos_validacion$grade_categ <- cut(datos_validacion$grade, breaks = c(0,4,9,13), labels = c(0,1,2))
datos_validacion$yr_renovated_catg <-cut(datos_validacion$yr_renovated, breaks=c(-0.5,1933, 2015), labels= c("0","1"))
datos_validacion$zipcode<-as.factor(datos_validacion$zipcode)
#codificar la variable Zipcode
datos_validacion$zona<-recode(datos_validacion$zipcode, " 98001=1; 98002=1; 98003=1; 98010=1; 98011=1; 98014=1; 98019=1; 98022=1; 98023=1; 98024=1; 98028=1; 98030=1; 98031=1; 98032=1; 98034=1; 98038=1; 98042=1; 98045=1; 98055=1; 98056=1; 98058=1; 98059=1; 98070=1; 98092=1 ;98106=1; 98108=1; 98118=1; 98125=1; 98126=1; 98133=1; 98144=1; 98146=1; 98148=1; 98155=1; 98166=1; 98168=1; 98178=1; 98188=1; 98198=1; 98004=2; 98005=2; 98006=2; 98007=2; 98008=2; 98027=2; 98029=2; 98033=2; 98039=2; 98040=2; 98052=2; 98053=2; 98065=2; 98072=2; 98074=2; 98075=2; 98077=2; 98102=2; 98103=2; 98105=2; 98107=2; 98109=2; 98112=2; 98115=2; 98116=2; 98117=2; 98119=2; 98122=2; 98136=2; 98177=2; 98199=2")
datos_validacion$zipcode = NULL
datos_validacion$bathrooms_group <- cut(datos_validacion$bathrooms_group, breaks = c(-1,2.5,8),labels=c(0,1))
datos_validacion_limpio <- datos_validacion[c(3,22:26,8:10,27,16,17,29,28,20,21)]
datos_validacion_limpio$log_above = NULL
datos_validacion_numeric <- datos_validacion_limpio %>% select_if(is.numeric)
#Aquí nos cargamos ya price_categ2
datos_train_limpio$price_categ2= NULL
datos_test_limpio$price_categ2= NULL
datos_validacion_limpio$price_categ2 = NULL
A continuación, se implementan dos métodos de agrupamiento no supervisado. Primero, un método jerárquico para averiguar (si es posible) la cantidad adecuada de clústeres y a continuación K-Means.
La agrupación es una técnica para agrupar puntos de datos similares en un grupo y separar las diferentes observaciones en diferentes grupos o grupos. En el Clustering Jerárquico los clusters se crean de manera que tengan un orden predeterminado. Existen En nuestro estudio se va a plicar un método aglomerativo que consiste en que cada observación se asigna a su propio clúster. Luego, se calcula la similitud (o distancia) entre cada uno de los clusters y los dos clusters más similares se fusionan en uno. Finalmente, los pasos 2 y 3 se repiten hasta que solo quede un grupo.
Aplicando este método a nuestros datos queremos observar si siguen algún patrón para poder agrupar las casas.
Escalamos las variables numéricas, es decir cada variable ahora tendrá una media cero y una desviación estándar uno.También se requiere los valores de distancia.La medida predeterminada para la función dist es ‘Euclidiana’.
datos_scale<- as.data.frame(scale(datos_train_numeric))
#Matriz de distancias
matriz_dist=dist(datos_scale)
En este caso particular, estimamos dos clústeres:
modelo_hc= hclust(matriz_dist, method = "average")
plot(modelo_hc)
rect.hclust(modelo_hc, k=2,border="red")
A continuación vemos como se han agrupado los datos marcando que el número de clústeres sea 2. Prácticamente todas las casas están en el grupo 1.
grupos2=cutree(modelo_hc,k=2)
table(datos_train_limpio$price_categ1, grupos2)
## grupos2
## 1 2
## B1 8724 10
## C1 6385 0
El método de K-Means basa su funcionamiento en agrupar los datos de entrada en un total de k conjuntos definidos por un centroide, cuya distancia con los puntos que pertenecen a cada uno de los datos es la menor posible.
Primero se van a realizar los gráficos para ver cómo están diferenciadas las casa por precio. Para poder visualizarlo en dos dimensiones se ha usado la función “prcomp” (PCA)
colores1= c("red","blue")
colores11 = colores1[datos_train_limpio$price_categ1]
plot(prcomp(datos_train_numeric, scale = T)$x[,1:2], type="n",main= "Dos categorías")
text(prcomp(datos_train_numeric, scale = T)$x[,1:2], as.character(datos_train_limpio$price_categ1), col=colores11)
A continuación se va a aplicar el método de agrupamiento de K-means, al igual que en el método anterior se le pasa la matriz de distacias y se van a agrupar los datos en dos conjuntos:
set.seed(1234)
modelkm1 <- kmeans(matriz_dist, centers=2)
#modelkm10 <- kmeans(matriz_dist, centers=10)
#table(modelkm10$cluster)
#modelkm2 <- kmeans(matriz_dist, centers=2)
#table(modelkm2$cluster) #asignación de observación a los cluster
# modelkm1$totss #Inercia total
# modelkm1$betweenss #Inercia inter grupos
# modelkm1$withinss #Inercia intra grupos
# modelkm1$tot.withinss #inercia total intra grupos
colores1= c("red","blue")
colores11 = colores1[datos_train_limpio$price_categ1]
plot(prcomp(datos_train_numeric, scale = T)$x[,1:2], type="n")
text(prcomp(datos_train_numeric, scale = T)$x[,1:2], as.character(modelkm1$cluster), col=colores11)
Se va a determinar la cantidad óptima de centroides a utilizar a partir del Método del Codo. Para ello, aplicaremos la función kmeans al conjunto de datos, variando en cada caso el valor de k, y acumulando los valores de WCSS obtenidos:
set.seed(1234)
wcss <- vector()
for(i in 1:20){
wcss[i] <- sum(kmeans(datos_scale, i)$withinss)
}
## Warning: did not converge in 10 iterations
ggplot() + geom_point(aes(x = 1:20, y = wcss), color = 'blue') +
geom_line(aes(x = 1:20, y = wcss), color = 'blue') +
ggtitle("Método del Codo") +
xlab('Cantidad de Centroides k') +
ylab('WCSS')
Como se observa en la gráfica, el K-óptimo que se podría aplicar sería de 10.
A continuación para visualizar si el agrupamiento que se ha llevado a cabo está relacionado con el precio de las casas(“Caras”, Baratas), se representa en el mapa que se muestra a continuación.
clusterkmeans=as.data.frame(modelkm1$cluster)
clusterkmeans$indice= as.integer(rownames(clusterkmeans))
colnames(clusterkmeans)[1]= "categoria_price_km"
clustering1= clusterkmeans[order(clusterkmeans$indice),]
Se visualiza como se reparten las casas según el agrupamiento que se ha realizado con K-Means:
center_lon = median(datos_train_limpio$long,na.rm = TRUE)
center_lat = median(datos_train_limpio$lat,na.rm = TRUE)
factpal1 <- colorFactor(c("green","red"),
clustering1$categoria_price_km )
leaflet(datos_train_limpio) %>% addProviderTiles("Esri.NatGeoWorldMap") %>%
addCircles(lng = ~long, lat = ~lat,
color = ~factpal1(clustering1$categoria_price_km)) %>%
# controls
setView(lng=center_lon, lat=center_lat,zoom = 12) %>%
addLegend("bottomright", pal = factpal1 , values = ~clustering1$categoria_price_km,
title = "Tipos de casas",
opacity = 1)
Cómo se observa en el mapa, y si lo comparamos con el de los datos iniciales, dista bastante. Por lo que deducimos que la agrupación que está relalizando K-Means no es muy buena.
K-medoids es un método de clustering muy similar a K-means en cuanto a que ambos agrupan las observaciones en K clusters, donde K es un valor preestablecido. La diferencia es que, en K-medoids, cada cluster está representado por una observación presente en el cluster (medoid),en nuestro estudio será una observación de una casa, mientras que en K-means cada cluster está representado por su centroide, que se corresponde con el promedio de todas las observaciones del cluster pero con ninguna en particular.
Este algoritmo es menos sensible al ruido y los valores atípicos, en comparación con k-means, porque usa medoides como centros de clúster en lugar de centroides. (utilizados en k-means).
# datos_train_limpio[,1:3]
# dist_eu <- as.matrix(dist(datos_train_limpio[,c(1,3)]))
# gower_daisy <- as.matrix(daisy(datos_train_limpio[,c(1,3)], metric = 'gower'))
# gower_daisy[1:10, 1:10]
datoskmedoids = datos_train_limpio[,-15]
medoids_model = pam(x = datoskmedoids, k = 2, keep.diss = TRUE, keep.data = TRUE)
medoids_model$medoids
## bedrooms bathrooms_group log_sqft_living log_lot sqft_basement_cat
## 9606 3 1 3.152288 3.870404 1
## 8735 4 2 3.406540 3.936111 2
## waterfront view condition grade_categ lat long zona
## 9606 1 1 3 2 47.4949 -122.239 1
## 8735 1 1 3 2 47.6477 -122.197 2
## yr_renovated_catg price_categ1
## 9606 1 1
## 8735 1 2
# plot(medoids_model, main = "silhouette plot")
grupos<-data.frame(datoskmedoids)
grupos<-cbind(grupos,data.frame(medoids_model$clustering))
grupos$medoids_model.clustering<-as.factor(grupos$medoids_model.clustering)
ggpairs(grupos,columns=c(1,2,3,15,12),mapping=aes(color=medoids_model.clustering))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
table(grupos$medoids_model.clustering)
##
## 1 2
## 8831 6288
#medoids_model$clustering
clustering= sort(medoids_model$clustering)
clustering=as.data.frame(medoids_model$clustering)
clustering$indice= as.integer(rownames(clustering))
colnames(clustering)[1]= "categoria_price"
clustering2= clustering[order(clustering$indice),]
center_lon = median(datoskmedoids$long,na.rm = TRUE)
center_lat = median(datoskmedoids$lat,na.rm = TRUE)
factpal2 <- colorFactor(c("green","red"),
clustering2$categoria_price )
leaflet(datoskmedoids) %>% addProviderTiles("Esri.NatGeoWorldMap") %>%
addCircles(lng = ~long, lat = ~lat,
color = ~factpal2(clustering2$categoria_price)) %>%
# controls
setView(lng=center_lon, lat=center_lat,zoom = 12) %>%
addLegend("bottomright", pal = factpal2 , values = ~clustering2$categoria_price,
title = "Tipos de casas",
opacity = 1)
Es un método que permite simplificar la complejidad de espacios muestrales con muchas dimensiones a la vez que conserva su información.
pca<-prcomp(datos_train_numeric,scale=T)
plot(prcomp(datos_train_numeric,scale=T))
#pca$center para saber la media de cada variable
#pca$scale para saber la desv.típica de cada variable
summary(prcomp(datos_train_numeric, scale=T))
## Importance of components:
## PC1 PC2 PC3 PC4 PC5
## Standard deviation 1.414 1.0885 0.9290 0.7850 0.57947
## Proportion of Variance 0.400 0.2369 0.1726 0.1232 0.06716
## Cumulative Proportion 0.400 0.6370 0.8096 0.9328 1.00000
biplot(x = pca, scale = 0, cex = 0.6, col = c("blue4", "brown3"))
Con este modelo se va a estudiar si existe relación entre el hecho de que una casa sea “cara” ó “barata” dependiendo de las características de las casas.Se va a generar un modelo en el que a apartir de las variables,…..prediga la probabilidad de que una casa sea barata o cara.
datos_train_rl <- datos_train_limpio[,-c(8,9)]
datos_train_rl$price_categ1<- recode(datos_train_rl$price_categ1, "'B1'=0; 'C1'=1")
#model_glm = glm(price ~., family = binomial, data =datos_train)
model_glm = glm(price_categ1 ~., family = binomial, data =datos_train_rl )
summary(model_glm)
##
## Call:
## glm(formula = price_categ1 ~ ., family = binomial, data = datos_train_rl)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.3536 -0.3998 -0.0883 0.3481 3.6096
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -759.91116 32.36457 -23.480 < 2e-16 ***
## bedrooms -0.23062 0.03834 -6.015 1.80e-09 ***
## bathrooms_group1 0.14717 0.06828 2.155 0.03113 *
## log_sqft_living 13.13156 0.32565 40.324 < 2e-16 ***
## log_lot 0.66619 0.08223 8.101 5.45e-16 ***
## sqft_basement_cat1 -0.45846 0.05867 -7.814 5.52e-15 ***
## waterfront1 1.35682 0.62801 2.161 0.03073 *
## view1 1.13634 0.22417 5.069 4.00e-07 ***
## view2 1.21640 0.14137 8.605 < 2e-16 ***
## view3 1.84563 0.21019 8.781 < 2e-16 ***
## view4 2.49482 0.49846 5.005 5.58e-07 ***
## lat 5.32574 0.24476 21.759 < 2e-16 ***
## long -3.75931 0.24366 -15.429 < 2e-16 ***
## zona2 3.33880 0.07036 47.451 < 2e-16 ***
## yr_renovated_catg1 0.37667 0.13792 2.731 0.00631 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 20592.9 on 15118 degrees of freedom
## Residual deviance: 8975.8 on 15104 degrees of freedom
## AIC: 9005.8
##
## Number of Fisher Scoring iterations: 6
predicciones <- ifelse(test = model_glm$fitted.values > 0.5, yes = 1, no = 0)
tabla_reglog <- table(model_glm$model$price_categ1, predicciones,
dnn = c("observaciones", "predicciones"))
#caras
pred.means2=tabla_reglog[2,2]/(tabla_reglog[2,2]+tabla_reglog[1,2])
rec.means2=tabla_reglog[2,2]/(tabla_reglog[2,2]+tabla_reglog[2,1])
F_caras_reglog=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_caras_reglog
## [1] 0.8477249
tabla_reglog
## predicciones
## observaciones 0 1
## 0 7813 921
## 1 1010 5375
#baratas
pred.means2=tabla_reglog[1,1]/(tabla_reglog[1,1]+tabla_reglog[2,1])
rec.means2=tabla_reglog[1,1]/(tabla_reglog[1,1]+tabla_reglog[1,2])
F_baratas_reglog=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_baratas_reglog
## [1] 0.8900154
tabla_reglog
## predicciones
## observaciones 0 1
## 0 7813 921
## 1 1010 5375
#F-Medida
F_reglog_train= (F_caras_reglog+F_baratas_reglog)/2
F_reglog_train
## [1] 0.8688702
accuracy_reglog_train = (tabla_reglog[1,1]+tabla_reglog[2,2]) / (tabla_reglog[1,1]+tabla_reglog[1,2]+tabla_reglog[2,1]+tabla_reglog[2,2])
accuracy_reglog_train
## [1] 0.8722799
Explicamos con palabras porqué este modelo no es bueno. En vez de seguir con test y validación se descartan.
#ggplot(subset(datos_train_rl, select = c(bedrooms,bathrooms_group, #log_sqft_living, log_lot, log_above, lat, long, price_categ1)), aes(x = #bathrooms_group,
#y = price_categ1)) + geom_point() + # geom_smooth(method = 'glm', method.args = list(family =
# 'binomial'), se=FALSE) +
#geom_line(data = data.frame(bathrooms_group = data.frame(bathrooms_group = #seq(1,8,1)), price_categ1 = predict(model_glm, data.frame(bathrooms_group = #seq(1,8,1)), type = "response")), colour = "navy")
Ver con cuantos vecinos:
suppressWarnings(suppressMessages(library(kknn)))
modelo <- train.kknn(price_categ1 ~ ., data = datos_train_limpio, kmax = 20)
modelo
##
## Call:
## train.kknn(formula = price_categ1 ~ ., data = datos_train_limpio, kmax = 20)
##
## Type of response variable: nominal
## Minimal misclassification: 0.1148886
## Best kernel: optimal
## Best k: 15
#k=1
#train
constante = 1
knn.train1=knn.cv(scale(datos_train_numeric),cl=as.factor(datos_train_limpio$price_categ1),k=constante)
tabla.knn.train1=table(knn.train1,datos_train_limpio$price_categ1)
pred.knn.train1=tabla.knn.train1[2,2]/(tabla.knn.train1[2,2]+tabla.knn.train1[1,2])
rec.knn.train1=tabla.knn.train1[2,2]/(tabla.knn.train1[2,2]+tabla.knn.train1[2,1])
F_medida.knn.train1=(5*pred.knn.train1*rec.knn.train1)/(4*pred.knn.train1+rec.knn.train1)
F_medida.knn.train1
## [1] 0.8367726
tabla.knn.train1
##
## knn.train1 B1 C1
## B1 7707 1073
## C1 1027 5312
#k=2
#train
knn_train2=knn.cv(scale(datos_train_numeric),cl=as.factor(datos_train_limpio$price_categ1),k=2)
tabla_knn_train2=table(knn_train2,datos_train_limpio$price_categ1)
pred_knn_train2=tabla_knn_train2[2,2]/(tabla_knn_train2[2,2]+tabla_knn_train2[1,2])
rec_knn_train2=tabla_knn_train2[2,2]/(tabla_knn_train2[2,2]+tabla_knn_train2[2,1])
F_medida_knn_train2=(5*pred_knn_train2*rec_knn_train2)/(4*pred_knn_train2+rec_knn_train2)
F_medida_knn_train2
## [1] 0.8304363
tabla_knn_train2
##
## knn_train2 B1 C1
## B1 7651 1082
## C1 1083 5303
center_lon = median(datos_train_limpio$long,na.rm = TRUE)
center_lat = median(datos_train_limpio$lat,na.rm = TRUE)
factpal1 <- colorFactor(c("green","red"),
knn_train2)
leaflet(datos_train_limpio) %>% addProviderTiles("Esri.NatGeoWorldMap") %>%
addCircles(lng = ~long, lat = ~lat,
color = ~factpal1(knn_train2)) %>%
# controls
setView(lng=center_lon, lat=center_lat,zoom = 12) %>%
addLegend("bottomright", pal = factpal1 , values = ~knn_train2,
title = "Precio (en miles de $)",
opacity = 1)
tabla_knn=table(datos_train_limpio$price_categ1,knn_train2)
#caras
pred_means_caras_knn=tabla_knn[2,2]/(tabla_knn[2,2]+tabla_knn[1,2])
rec_means_knn=tabla_knn[2,2]/(tabla_knn[2,2]+tabla_knn[2,1])
F_caras_knn_train=(2*pred_means_caras_knn*rec_means_knn)/(pred_means_caras_knn+rec_means_knn)
F_caras_knn_train
## [1] 0.8304753
tabla_knn
## knn_train2
## B1 C1
## B1 7651 1083
## C1 1082 5303
#baratas
pred_means_baratas_knn=tabla_knn[1,1]/(tabla_knn[1,1]+tabla_knn[2,1])
rec_means_baratas_knn=tabla_knn[1,1]/(tabla_knn[1,1]+tabla_knn[1,2])
F_baratas_knn_train=(2*pred_means_baratas_knn*rec_means_baratas_knn)/(pred_means_baratas_knn+rec_means_baratas_knn)
F_baratas_knn_train
## [1] 0.876052
tabla_knn
## knn_train2
## B1 C1
## B1 7651 1083
## C1 1082 5303
#Media de la F-MEDIDA para KNN
F_knn_train= (F_caras_knn_train+F_baratas_knn_train)/2
F_knn_train
## [1] 0.8532636
accuracy_knn_train= (tabla_knn[1,1]+tabla_knn[2,2]) / (tabla_knn[1,1]+tabla_knn[1,2]+tabla_knn[2,1]+tabla_knn[2,2])
accuracy_knn_train
## [1] 0.8568027
knn_test2=knn(scale(datos_train_numeric),scale(datos_test_numeric),cl=datos_train_limpio$price_categ1,k=2)
tabla2_test=table(datos_test_limpio$price_categ1,knn_test2)
#caras
pred.means2=tabla2_test[2,2]/(tabla2_test[2,2]+tabla2_test[1,2])
rec.means2=tabla2_test[2,2]/(tabla2_test[2,2]+tabla2_test[2,1])
F_medida.means2=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_medida.means2
## [1] 0.8109322
tabla2_test
## knn_test2
## B1 C1
## B1 1653 258
## C1 247 1083
#baratas
pred.means2=tabla2_test[1,1]/(tabla2_test[1,1]+tabla2_test[2,1])
rec.means2=tabla2_test[1,1]/(tabla2_test[1,1]+tabla2_test[1,2])
F_medida.means2=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_medida.means2
## [1] 0.8674888
tabla2_test
## knn_test2
## B1 C1
## B1 1653 258
## C1 247 1083
accuracy = (tabla2_test[1,1]+tabla2_test[2,2]) / (tabla2_test[1,1]+tabla2_test[1,2]+tabla2_test[2,1]+tabla2_test[2,2])
accuracy
## [1] 0.8441839
knn.val=knn(scale(datos_train_numeric),scale(datos_validacion_numeric),cl=datos_train_limpio$price_categ1,k=2)
tabla2_val=table(datos_validacion_limpio$price_categ1,knn.val)
#caras
pred.means2=tabla2_val[2,2]/(tabla2_val[2,2]+tabla2_val[1,2])
rec.means2=tabla2_val[2,2]/(tabla2_val[2,2]+tabla2_val[2,1])
F_medida.means2=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_medida.means2
## [1] 0.8139274
tabla2_val
## knn.val
## B1 C1
## B1 1659 247
## C1 250 1087
#baratas
pred.means2=tabla2_val[1,1]/(tabla2_val[1,1]+tabla2_val[2,1])
rec.means2=tabla2_val[1,1]/(tabla2_val[1,1]+tabla2_val[1,2])
F_medida.means2=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_medida.means2
## [1] 0.8697248
tabla2_val
## knn.val
## B1 C1
## B1 1659 247
## C1 250 1087
accuracy = (tabla2_val[1,1]+tabla2_val[2,2]) / (tabla2_val[1,1]+tabla2_val[1,2]+tabla2_val[2,1]+tabla2_val[2,2])
accuracy
## [1] 0.8467468
#quitamos lat y long:
datos_decision_tree<-datos_train_limpio[,-c(11,12)]
decisiontree_model=rpart(price_categ1~., data=datos_decision_tree,parms=list(split="gini") )
#resumen=summary(decisiontree_model)
print(decisiontree_model)
## n= 15119
##
## node), split, n, loss, yval, (yprob)
## * denotes terminal node
##
## 1) root 15119 6385 B1 (0.57768371 0.42231629)
## 2) log_sqft_living< 3.35698 9914 2474 B1 (0.75045390 0.24954610)
## 4) lat< 47.55595 4723 252 B1 (0.94664408 0.05335592) *
## 5) lat>=47.55595 5191 2222 B1 (0.57195145 0.42804855)
## 10) lat>=47.69045 2056 353 B1 (0.82830739 0.17169261) *
## 11) lat< 47.69045 3135 1266 C1 (0.40382775 0.59617225)
## 22) log_sqft_living< 3.173768 1313 467 B1 (0.64432597 0.35567403) *
## 23) log_sqft_living>=3.173768 1822 420 C1 (0.23051592 0.76948408) *
## 3) log_sqft_living>=3.35698 5205 1294 C1 (0.24860711 0.75139289)
## 6) lat< 47.51495 1422 451 B1 (0.68284107 0.31715893)
## 12) log_lot< 4.175105 973 183 B1 (0.81192189 0.18807811) *
## 13) log_lot>=4.175105 449 181 C1 (0.40311804 0.59688196) *
## 7) lat>=47.51495 3783 323 C1 (0.08538197 0.91461803) *
#kk=as.data.frame(summary(decisiontree_model))
fancyRpartPlot(decisiontree_model)
#rpart.plot(decisiontree_model)
#prp(decisiontree_model, type = 2, extra = 6, split.font = 3)
tabla_train_arbol=table(obs = datos_decision_tree$price_categ1, pred = predict(decisiontree_model, type = "class"))
tabla_train_arbol
## pred
## obs B1 C1
## B1 7810 924
## C1 1255 5130
#caras
pred.means2=tabla_train_arbol[2,2]/(tabla_train_arbol[2,2]+tabla_train_arbol[1,2])
rec.means2=tabla_train_arbol[2,2]/(tabla_train_arbol[2,2]+tabla_train_arbol[2,1])
F_caras_dt_train=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_caras_dt_train
## [1] 0.8248251
tabla_train_arbol
## pred
## obs B1 C1
## B1 7810 924
## C1 1255 5130
#baratas
pred.means2=tabla_train_arbol[1,1]/(tabla_train_arbol[1,1]+tabla_train_arbol[2,1])
rec.means2=tabla_train_arbol[1,1]/(tabla_train_arbol[1,1]+tabla_train_arbol[1,2])
F_baratas_dt_train=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_baratas_dt_train
## [1] 0.8775774
tabla_train_arbol
## pred
## obs B1 C1
## B1 7810 924
## C1 1255 5130
#Media de la F-MEDIDA para Decision Tree
F_dt_train= (F_caras_dt_train+F_baratas_dt_train)/2
F_dt_train
## [1] 0.8512013
accuracy_dt_train = (tabla_train_arbol[1,1]+tabla_train_arbol[2,2]) / (tabla_train_arbol[1,1]+tabla_train_arbol[1,2]+tabla_train_arbol[2,1]+tabla_train_arbol[2,2])
accuracy_dt_train
## [1] 0.8558767
tabla.test.arbol1=table(obs = datos_test_limpio$price_categ1, pred = predict(decisiontree_model,datos_test_limpio[,-15], type ="class"))
tabla.test.arbol1
## pred
## obs B1 C1
## B1 1690 221
## C1 262 1068
#caras
pred.means2=tabla.test.arbol1[2,2]/(tabla.test.arbol1[2,2]+tabla.test.arbol1[1,2])
rec.means2=tabla.test.arbol1[2,2]/(tabla.test.arbol1[2,2]+tabla.test.arbol1[2,1])
F_medida.means2=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_medida.means2
## [1] 0.8155785
tabla.test.arbol1
## pred
## obs B1 C1
## B1 1690 221
## C1 262 1068
#baratas
pred.means2=tabla.test.arbol1[1,1]/(tabla.test.arbol1[1,1]+tabla.test.arbol1[2,1])
rec.means2=tabla.test.arbol1[1,1]/(tabla.test.arbol1[1,1]+tabla.test.arbol1[1,2])
F_medida.means2=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_medida.means2
## [1] 0.8749676
tabla.test.arbol1
## pred
## obs B1 C1
## B1 1690 221
## C1 262 1068
accuracy = (tabla.test.arbol1[1,1]+tabla.test.arbol1[2,2]) / (tabla.test.arbol1[1,1]+tabla.test.arbol1[1,2]+tabla.test.arbol1[2,1]+tabla.test.arbol1[2,2])
accuracy
## [1] 0.8509719
tabla.validacion.arbol1=table(obs = datos_validacion_limpio$price_categ1, pred = predict(decisiontree_model,datos_validacion_limpio[,-15], type ="class"))
tabla.validacion.arbol1
## pred
## obs B1 C1
## B1 1700 206
## C1 263 1074
#caras
pred.means2=tabla.validacion.arbol1[2,2]/(tabla.validacion.arbol1[2,2]+tabla.validacion.arbol1[1,2])
rec.means2=tabla.validacion.arbol1[2,2]/(tabla.validacion.arbol1[2,2]+tabla.validacion.arbol1[2,1])
F_medida.means2=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_medida.means2
## [1] 0.8207872
tabla.validacion.arbol1
## pred
## obs B1 C1
## B1 1700 206
## C1 263 1074
#baratas
pred.means2=tabla.validacion.arbol1[1,1]/(tabla.validacion.arbol1[1,1]+tabla.validacion.arbol1[2,1])
rec.means2=tabla.validacion.arbol1[1,1]/(tabla.validacion.arbol1[1,1]+tabla.validacion.arbol1[1,2])
F_medida.means2=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_medida.means2
## [1] 0.87878
tabla.validacion.arbol1
## pred
## obs B1 C1
## B1 1700 206
## C1 263 1074
accuracy = (tabla.validacion.arbol1[1,1]+tabla.validacion.arbol1[2,2]) / (tabla.validacion.arbol1[1,1]+tabla.validacion.arbol1[1,2]+tabla.validacion.arbol1[2,1]+tabla.validacion.arbol1[2,2])
accuracy
## [1] 0.8553808
randomforest_model=randomForest(price_categ1~., data=datos_train_limpio,parms=list(split="gini") )
#print(randomforest_model)
tabla_randomforest=table(obs = datos_train_limpio$price_categ1, pred = predict(randomforest_model, type = "class") )
tabla_randomforest
## pred
## obs B1 C1
## B1 8001 733
## C1 726 5659
#caras
pred.means2=tabla_randomforest[2,2]/(tabla_randomforest[2,2]+tabla_randomforest[1,2])
rec.means2=tabla_randomforest[2,2]/(tabla_randomforest[2,2]+tabla_randomforest[2,1])
F_caras_rf_train=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_caras_rf_train
## [1] 0.8858104
tabla_randomforest
## pred
## obs B1 C1
## B1 8001 733
## C1 726 5659
#baratas
pred.means2=tabla_randomforest[1,1]/(tabla_randomforest[1,1]+tabla_randomforest[2,1])
rec.means2=tabla_randomforest[1,1]/(tabla_randomforest[1,1]+tabla_randomforest[1,2])
F_baratas_rf_train=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_baratas_rf_train
## [1] 0.9164424
tabla_randomforest
## pred
## obs B1 C1
## B1 8001 733
## C1 726 5659
#Media de la F-MEDIDA para Random Forest
F_rf_train= (F_caras_rf_train+F_baratas_rf_train)/2
F_rf_train
## [1] 0.9011264
accuracy_rf_train = (tabla_randomforest[1,1]+tabla_randomforest[2,2]) / (tabla_randomforest[1,1]+tabla_randomforest[1,2]+tabla_randomforest[2,1]+tabla_randomforest[2,2])
accuracy_rf_train
## [1] 0.9034989
tabla_randomforest_test=table(obs = datos_test_limpio$price_categ1, pred = predict(randomforest_model,datos_test_limpio[,-15], type ="class"))
tabla_randomforest_test
## pred
## obs B1 C1
## B1 1745 166
## C1 144 1186
#caras
pred.means2=tabla_randomforest_test[2,2]/(tabla_randomforest_test[2,2]+tabla_randomforest_test[1,2])
rec.means2=tabla_randomforest_test[2,2]/(tabla_randomforest_test[2,2]+tabla_randomforest_test[2,1])
F_medida.means2=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_medida.means2
## [1] 0.8844146
tabla_randomforest_test
## pred
## obs B1 C1
## B1 1745 166
## C1 144 1186
#baratas
pred.means2=tabla_randomforest_test[1,1]/(tabla_randomforest_test[1,1]+tabla_randomforest_test[2,1])
rec.means2=tabla_randomforest_test[1,1]/(tabla_randomforest_test[1,1]+tabla_randomforest_test[1,2])
F_medida.means2=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_medida.means2
## [1] 0.9184211
tabla_randomforest_test
## pred
## obs B1 C1
## B1 1745 166
## C1 144 1186
accuracy = (tabla_randomforest_test[1,1]+tabla_randomforest_test[2,2]) / (tabla_randomforest_test[1,1]+tabla_randomforest_test[1,2]+tabla_randomforest_test[2,1]+tabla_randomforest_test[2,2])
accuracy
## [1] 0.9043505
tabla_randomforest_val=table(obs = datos_validacion_limpio$price_categ1, pred = predict(randomforest_model,datos_validacion_limpio[,-15], type ="class"))
tabla_randomforest_val
## pred
## obs B1 C1
## B1 1756 150
## C1 148 1189
#caras
pred.means2=tabla_randomforest_val[2,2]/(tabla_randomforest_val[2,2]+tabla_randomforest_val[1,2])
rec.means2=tabla_randomforest_val[2,2]/(tabla_randomforest_val[2,2]+tabla_randomforest_val[2,1])
F_medida.means2=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_medida.means2
## [1] 0.8886398
tabla_randomforest_val
## pred
## obs B1 C1
## B1 1756 150
## C1 148 1189
#baratas
pred.means2=tabla_randomforest_val[1,1]/(tabla_randomforest_val[1,1]+tabla_randomforest_val[2,1])
rec.means2=tabla_randomforest_val[1,1]/(tabla_randomforest_val[1,1]+tabla_randomforest_val[1,2])
F_medida.means2=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_medida.means2
## [1] 0.9217848
tabla_randomforest_val
## pred
## obs B1 C1
## B1 1756 150
## C1 148 1189
accuracy = (tabla_randomforest_val[1,1]+tabla_randomforest_val[2,2]) / (tabla_randomforest_val[1,1]+tabla_randomforest_val[1,2]+tabla_randomforest_val[2,1]+tabla_randomforest_val[2,2])
accuracy
## [1] 0.9081098
modelo_svm <- e1071::svm(formula = price_categ1 ~., data = datos_train_limpio, kernel = "linear" , probability =TRUE)
summary(modelo_svm)
##
## Call:
## svm(formula = price_categ1 ~ ., data = datos_train_limpio, kernel = "linear",
## probability = TRUE)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: linear
## cost: 1
##
## Number of Support Vectors: 4702
##
## ( 2352 2350 )
##
##
## Number of Classes: 2
##
## Levels:
## B1 C1
#plot(modelo_svm, datos_train_limpio_dt)
#set.seed(1)
#tune.out = tune(svm,price_categ1 ~ ., data = datos_train_limpio_dt, kernel = "linear", ranges = list(cost = c(0.001,
#0.01, 0.1, 1, 5, 10, 100)))
tabla_svm_train=table(obs = datos_train_limpio$price_categ1, pred = predict(modelo_svm,datos_train_limpio[,-15], type ="class"))
tabla_svm_train
## pred
## obs B1 C1
## B1 7835 899
## C1 1034 5351
#caras
pred.means2=tabla_svm_train[2,2]/(tabla_svm_train[2,2]+tabla_svm_train[1,2])
rec.means2=tabla_svm_train[2,2]/(tabla_svm_train[2,2]+tabla_svm_train[2,1])
F_caras_svm_train=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_caras_svm_train
## [1] 0.8470123
tabla_svm_train
## pred
## obs B1 C1
## B1 7835 899
## C1 1034 5351
#baratas
pred_means2=tabla_svm_train[1,1]/(tabla_svm_train[1,1]+tabla_svm_train[2,1])
rec_means2=tabla_svm_train[1,1]/(tabla_svm_train[1,1]+tabla_svm_train[1,2])
F_baratas_svm_train=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_baratas_svm_train
## [1] 0.8470123
tabla_svm_train
## pred
## obs B1 C1
## B1 7835 899
## C1 1034 5351
#F-Medida para SVM
F_svm_train= (F_caras_svm_train+F_baratas_svm_train)/2
F_svm_train
## [1] 0.8470123
accuracy_svm_train = (tabla_svm_train[1,1]+tabla_svm_train[2,2]) / (tabla_svm_train[1,1]+tabla_svm_train[1,2]+tabla_svm_train[2,1]+tabla_svm_train[2,2])
accuracy_svm_train
## [1] 0.8721476
tabla_svm_test=table(obs = datos_test_limpio$price_categ1, pred = predict(modelo_svm,datos_test_limpio[,-15], type ="class"))
tabla_svm_test
## pred
## obs B1 C1
## B1 1694 217
## C1 213 1117
#caras
pred.means2=tabla_svm_test[2,2]/(tabla_svm_test[2,2]+tabla_svm_test[1,2])
rec.means2=tabla_svm_test[2,2]/(tabla_svm_test[2,2]+tabla_svm_test[2,1])
F_medida.means2=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_medida.means2
## [1] 0.8385886
tabla_svm_test
## pred
## obs B1 C1
## B1 1694 217
## C1 213 1117
#baratas
pred.means2=tabla_svm_test[1,1]/(tabla_svm_test[1,1]+tabla_svm_test[2,1])
rec.means2=tabla_svm_test[1,1]/(tabla_svm_test[1,1]+tabla_svm_test[1,2])
F_medida.means2=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_medida.means2
## [1] 0.8873756
tabla_svm_test
## pred
## obs B1 C1
## B1 1694 217
## C1 213 1117
accuracy = (tabla_svm_test[1,1]+tabla_svm_test[2,2]) / (tabla_svm_test[1,1]+tabla_svm_test[1,2]+tabla_svm_test[2,1]+tabla_svm_test[2,2])
accuracy
## [1] 0.8673249
tabla_svm_val=table(obs = datos_validacion_limpio$price_categ1, pred = predict(modelo_svm,datos_validacion_limpio[,-15], type ="class"))
tabla_svm_val
## pred
## obs B1 C1
## B1 1711 195
## C1 206 1131
#caras
pred.means2=tabla_svm_val[2,2]/(tabla_svm_val[2,2]+tabla_svm_val[1,2])
rec.means2=tabla_svm_val[2,2]/(tabla_svm_val[2,2]+tabla_svm_val[2,1])
F_medida.means2=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_medida.means2
## [1] 0.8494179
#baratas
pred.means2=tabla_svm_val[1,1]/(tabla_svm_val[1,1]+tabla_svm_val[2,1])
rec.means2=tabla_svm_val[1,1]/(tabla_svm_val[1,1]+tabla_svm_val[1,2])
F_medida.means2=(2*pred.means2*rec.means2)/(pred.means2+rec.means2)
F_medida.means2
## [1] 0.8951086
tabla_svm_val
## pred
## obs B1 C1
## B1 1711 195
## C1 206 1131
accuracy = (tabla_svm_val[1,1]+tabla_svm_val[2,2]) / (tabla_svm_val[1,1]+tabla_svm_val[1,2]+tabla_svm_val[2,1]+tabla_svm_val[2,2])
accuracy
## [1] 0.8763491
datos_train_gam <- datos_train_limpio
datos_train_gam$price <- datos_train$price
model_gam <- gam(price ~ s(bedrooms, bs = "cr") + bathrooms_group + s(log_sqft_living, bs = "cr") +
s(log_lot, bs = "cr") + sqft_basement_cat + waterfront + view + s(lat, bs = "cr") +
s(long, bs = "cr") + zona + yr_renovated_catg, data = datos_train_gam)
summary(model_gam)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## price ~ s(bedrooms, bs = "cr") + bathrooms_group + s(log_sqft_living,
## bs = "cr") + s(log_lot, bs = "cr") + sqft_basement_cat +
## waterfront + view + s(lat, bs = "cr") + s(long, bs = "cr") +
## zona + yr_renovated_catg
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 460094 3772 121.978 < 2e-16 ***
## bathrooms_group1 35433 4091 8.662 < 2e-16 ***
## sqft_basement_cat1 -33249 3315 -10.029 < 2e-16 ***
## waterfront1 523178 20514 25.503 < 2e-16 ***
## view1 90905 11827 7.686 1.61e-14 ***
## view2 82058 7077 11.595 < 2e-16 ***
## view3 155689 9712 16.031 < 2e-16 ***
## view4 287662 14672 19.606 < 2e-16 ***
## zona2 126977 5000 25.398 < 2e-16 ***
## yr_renovated_catg1 58719 7070 8.306 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(bedrooms) 8.739 8.957 19.8 <2e-16 ***
## s(log_sqft_living) 8.911 8.997 1509.6 <2e-16 ***
## s(log_lot) 8.339 8.802 37.3 <2e-16 ***
## s(lat) 8.872 8.995 280.6 <2e-16 ***
## s(long) 8.935 8.999 171.2 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.781 Deviance explained = 78.2%
## GCV = 3.0207e+10 Scale est. = 3.0099e+10 n = 15119
plot(model_gam, ylab = "price")
par(mfrow=c(3,3))
visreg(model_gam)
datos_test_limpio$price <- datos_test$price
datos_test_limpio$price_predict = predict(model_gam, datos_test_limpio)
datos_test_limpio$precio_diff <- abs(datos_test_limpio$price - datos_test_limpio$price_predict)
dif_media_test <- mean(datos_test_limpio$precio_diff)
dif_media_test
## [1] 108157.9
datos_validacion_limpio$price <- datos_validacion$price
datos_validacion_limpio$price_predict = predict(model_gam, datos_validacion_limpio)
datos_validacion_limpio$precio_diff <- abs(datos_validacion_limpio$price - datos_validacion_limpio$price_predict)
dif_media_test <- mean(datos_validacion_limpio$precio_diff)
dif_media_test
## [1] 105924.8
#knn.train2=knn.cv(scale(datos_train_numeric),cl=as.factor(datos_train_limpio$pri#ce_categ1),k=2)
# x <- iris[,-5]
#y <- iris[,5]
obj2 <- tune.knn(x=scale(datos_train_numeric),
y=as.factor(datos_train_limpio$price_categ1), k = 1:10,
tunecontrol = tune.control(sampling = "boot"))
#summary(obj2)
plot(obj2)
#n^(4/(d+4))
#n es el número de datos de entrenamiento
#d es la dimensión de los datos
n=dim(datos_train_numeric)[1]
d=ncol(datos_train_numeric)
k_optimo = as.integer(n^(4/(d+4)))
k_optimo
## [1] 72
rango=1:(k_optimo+10)
modelos=c()
f1_modelos=c()
for (i in rango){
model_knn=knn.cv(scale(datos_train_numeric),cl=as.factor(datos_train_limpio$price_categ1),k=i)
tabla=table(datos_train_limpio$price_categ1,model_knn)
# f1 casas caras
pred_means_caras=tabla[2,2]/(tabla[2,2]+tabla[1,2])
rec_means_caras=tabla[2,2]/(tabla[2,2]+tabla[2,1])
f1_caras=(2*pred_means_caras*rec_means_caras)/(pred_means_caras+rec_means_caras)
# f1 casas baratas
pred_means_baratas=tabla[1,1]/(tabla[1,1]+tabla[2,1])
rec_means_baratas=tabla[1,1]/(tabla[1,1]+tabla[1,2])
f1_baratas=(2*pred_means_baratas*rec_means_baratas)/(pred_means_baratas+rec_means_baratas)
f1_total = (f1_baratas + f1_caras)/2
f1_modelos=c(f1_modelos,f1_total)
}
f1_modelos
## [1] 0.8575267 0.8503791 0.8735162 0.8683478 0.8753198 0.8755879 0.8766022
## [8] 0.8751717 0.8786105 0.8752529 0.8780367 0.8790804 0.8792161 0.8771628
## [15] 0.8787699 0.8775224 0.8797486 0.8789995 0.8782733 0.8796625 0.8784356
## [22] 0.8782627 0.8793153 0.8790047 0.8790228 0.8788375 0.8785713 0.8780670
## [29] 0.8787176 0.8775505 0.8775171 0.8767733 0.8767294 0.8764420 0.8769355
## [36] 0.8766922 0.8764966 0.8755295 0.8766043 0.8753725 0.8760616 0.8751011
## [43] 0.8745799 0.8743301 0.8743193 0.8731636 0.8740480 0.8743576 0.8743301
## [50] 0.8736573 0.8742159 0.8732341 0.8733860 0.8729248 0.8728759 0.8729464
## [57] 0.8729518 0.8716501 0.8719591 0.8709986 0.8711076 0.8710807 0.8706303
## [64] 0.8711186 0.8708801 0.8710535 0.8701859 0.8700173 0.8692036 0.8699356
## [71] 0.8698332 0.8701804 0.8694319 0.8698662 0.8689657 0.8694485 0.8691829
## [78] 0.8696820 0.8687928 0.8678551 0.8685051 0.8676693
rango
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
## [51] 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
## [76] 76 77 78 79 80 81 82
pos_max <- which.max(f1_modelos)
rango[pos_max]
## [1] 17
Una vez que se han entrenado y optimizado distintos modelos, se tiene que identificar cuál de ellos consigue mejores resultados para el problema en cuestión, en este caso, predecir si una casa es barata o cara. Con los datos disponibles, existen dos formas de comparar los modelos. Si bien las dos no tienen por qué dar los mismos resultados, son complementarias a la hora de tomar una decisión final.
models_cross = data.frame(
"modelo"= c('GLM','knn','Decision_Tree','Random_Forest','SVM'),
"F_Medida_train" = c(F_reglog_train,F_knn_train,F_dt_train,F_rf_train,F_svm_train))
plot(x = models_cross$modelo, y= models_cross$F_Medida_train,fill=models_cross$modelo)
ggplot(data=models_cross, aes(x=modelo, y=F_Medida_train, fill=modelo)) +
geom_bar(stat="identity", position="dodge")
# REGRESIÓN LOGÍSTICA
predictions_glm <- predict(model_glm, newdata = datos_train_rl, type = "response")
pred_glm <- prediction(predictions_glm, datos_train_rl$price_categ1)
perf_glm <- performance(pred_glm,"tpr","fpr")
# ARBOL DE DECISION
predictions_tree <- predict(decisiontree_model, newdata = datos_decision_tree, type = "prob")
pred_tree = prediction(predictions_tree[,2], datos_decision_tree$price_categ1)
pref_tree = performance(pred_tree, "tpr", "fpr")
# RANDOM FOREST
predictions_rf <- predict(randomforest_model, new_data=datos_train_limpio, type = "prob")
pred_rf <- prediction(predictions_rf[,2],datos_train_limpio$price_categ1)
perf_rf <- performance(pred_rf,"tpr","fpr")
# SVM
predictions_svm <- predict(modelo_svm, newdata=datos_train_limpio, probability = TRUE)
prob_svm<-attr(predictions_svm,"probabilities")
pred_svm <- prediction(prob_svm[,2], datos_train_limpio$price_categ)
perf_svm <- performance(pred_svm,"tpr","fpr")
plot(perf_glm,col="blue")
plot(pref_tree, col="red",add = TRUE)
plot(perf_rf,col="green", add = TRUE)
plot(perf_svm, col="yellow",add = TRUE)
legend(x="right" ,legend=c("GLM","DT","RF","SVM"),
fill=c("blue","red","green","yellow"),cex=0.8)
cambiar el modelo para que aumente 1000 más la variable respuesta
Primero transformaciones de train, test validacion Luego cada modelo por separado con su train y test y validacion